function [alphas,betas,R2s, mats_yrs] = CSregs_real(gesol, Phi_Nz, stationary_prob_Nz)
% INPUTS:
% 1. [nominalylds_mean, nominalylds_volatilities, nominalbondprice_components] = nominal_logyld_moments_accelerate( F, parms, gesol, ...
%    inflation_data_est, test_mats, Ygrid, stationary_prob_NzY);
% 2. [Ygrid, Phi_NzY] = TransitionMatrix_NzY_Sparse(parms, gesol, rho_pi, F, dYmin);
% 3. stationary_prob_NzY
% OUTPUTS:
% campbell-shiller regressions for n = 2, 3, 4, 5 years
% Auxiliary functions used:
% C = calc_cov_NzYXe(Ygrid,Phi_NzY,stationary_prob_NzY,inflation_data_est,...
%    F1_Nz,F1_Y,F1_X,F1_eps,F2_Nz,F2_Y,F2_X,F2_eps,h)
% V = calc_var_NzYXe(Ygrid,stationary_prob_NzY,inflation_data_est,...
%     F_Nz,F_Y,F_X,F_eps)
    
    mats_yrs = 2:5;
    
    alphas = zeros(size(mats_yrs));
    betas = zeros(size(mats_yrs));
    R2s = zeros(size(mats_yrs));
    
    for i = 1:numel(mats_yrs)
        
        mat_yr = mats_yrs(i);
        mat_mth = mat_yr*12;
        
        idx_1y = find(gesol.real_bonds.maturities==12);
        idx_ny = find(gesol.real_bonds.maturities==12*mat_yr);
        idx_nless1y = find(gesol.real_bonds.maturities==12*(mat_yr-1));
        
        if isempty(idx_1y) || isempty(idx_ny) || isempty(idx_nless1y)
            error('Some yields are not available');
        end
        
        % compute coefficients for RHS: slope = (y(n) - y(1))/(n-1)
        rhs_Nz = (-log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny)))/mat_yr ...
            + log(squeeze(gesol.real_bonds.prices(:,:,:,idx_1y))))/(mat_yr - 1);
        
        % compute coefficients for LHS = y(n-1,t+1) - y(n,t)
        % lhs1 = -y(n,t), lhs2 at h=12 is y(n-1,t+1)
        
        lhs1_Nz = log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny)))/mat_yr;
        lhs2_Nz = -log(squeeze(gesol.real_bonds.prices(:,:,:,idx_nless1y)))/(mat_yr-1);
        
        % reg coeffs
        mean_rhs = sum(stationary_prob_Nz(:).*rhs_Nz(:));
        mean_lhs1 = sum(stationary_prob_Nz(:).*lhs1_Nz(:));
        mean_lhs2 = sum(stationary_prob_Nz(:).*lhs2_Nz(:));
        var_rhs = sum(stationary_prob_Nz(:).*(rhs_Nz(:).^2)) - mean_rhs^2;
        var_lhs1 = sum(stationary_prob_Nz(:).*(lhs1_Nz(:).^2)) - mean_lhs1^2;
        var_lhs2 = sum(stationary_prob_Nz(:).*(lhs2_Nz(:).^2)) - mean_lhs2^2;
        cov_lhs1_lhs2 = sum(stationary_prob_Nz(:).*lhs1_Nz(:).*((Phi_Nz^12)*lhs2_Nz(:))) - mean_lhs1*mean_lhs2;
        var_lhs = var_lhs1 + var_lhs2 + 2*cov_lhs1_lhs2;
        cov_rhs_lhs1 = sum(stationary_prob_Nz(:).*rhs_Nz(:).*lhs1_Nz(:)) - mean_rhs*mean_lhs1;
        cov_rhs_lhs2 = sum(stationary_prob_Nz(:).*rhs_Nz(:).*((Phi_Nz^12)*lhs2_Nz(:))) - mean_rhs*mean_lhs2;
        
        betas(i) = (cov_rhs_lhs1 + cov_rhs_lhs2)/var_rhs;
        alphas(i) = mean_lhs1 + mean_lhs2 - betas(i)*mean_rhs;
        R2s(i) = (betas(i)^2)*var_rhs/var_lhs;
        
    end

end